Task 1.6

Find task description and steps here

Find and download the source code from here

The model (metarange) predicts abundance of the species for all the cells in the landscape (in this case two states in Germany - Bavaria and Baden-Württemberg).

Step 1 - Import and merge all output files

all_data <- list.files(recursive = T) %>% #get a list of the paths to datasets files
  map_dfr(function(path) {
 #retrieve any characters from the beginning of the path until /
  foldername <- str_extract(path, "^.+(?=/)")
  fread(path) %>% #import the files  
  mutate(folder = foldername) # add a new collumn with the folder name
})
invisible(gc())

all_data <- all_data %>%
  separate(folder, c("date", "time", "species", "scenario"), sep = "_", remove = FALSE) %>% #takes a VERY long time ^(more than 40 minutes)
  dplyr::select(-c("folder", "date", "time"))
invisible(gc())

#check for negative values in abundance
summary(all_data)
#rows with negative abundances
bb <- as.data.frame(which(all_data < 0 , arr.ind=TRUE))
bb <- bb[bb$col == 4,]

#values <- as.vector(bb$row)
aa <- all_data[values,]

#write tsv file with all species in all scenarios
write_tsv(all_data, "full_dataset.tsv")

#Import the complete dataset 
all_data <- fread("./full_dataset.tsv", integer64 = "numeric")
nvisible(gc())

Step 2 - Create new variables

Occupancy - cells with at least one individual

The following variables are all calculated in relation to time step 25, which corresponds to the end of the burn-in period that allows the species to stabilize (the model has reached it stationary region).

Abundance change - \(abundance_{tn} - abundance_{t25}\)

Occupancy change - \(ocupancy_{tn} - occupancy_{t25}\)

Habitat change - \(habitat\ suitability_{tn} - habitat\ suitability_{t25}\)

#################################
# STEP 2 - CREATE NEW VARIABLES #
#################################

# occupancy
all_data <- all_data %>%
  mutate(occupancy = ifelse(all_data$abundance > 0, 1, 0)) 
invisible(gc())

# abundance change
all_data <- all_data %>%
  group_by(species, scenario) %>% 
  mutate(abund_change <- abundance - abundance[t == 25]) %>%
  ungroup()
invisible(gc())

# occupancy change
all_data <- all_data %>%
  group_by(species, scenario) %>% 
  mutate(occup_change <- occupancy - occupancy[t == 25]) %>%
  ungroup()
invisible(gc())

# habitat suitability change
all_data <- all_data %>%
  group_by(species, scenario) %>% 
  mutate(habitat_change <- habitat - habitat[t == 25]) %>%
  ungroup()
invisible(gc())

Step 3 - Population level plots (per species)

Population-level = value for the entire population (e.g. average of cell values for the entire population in a given time step)

For the following plots the burn in period is not included in the graphical representation, since it corresponds to the preliminary steps of the model.

##########################################
# STEP 3 - CREATE POPULATION LEVEL PLOTS #
##########################################

# calculate mean values
data_all_mean_stdev<-all_data[, c(1,4:10)]

# mean and standard deviation for ABUNDANCE, REPRODUCTION, CARRYING CAPACITY, HABITAT SUITABILITY and MORTALITY variables
mean_stDEV <- data_all_mean_stdev %>%
  dplyr::group_by(species, scenario, t) %>%
  dplyr::summarise(across(c(abundance, reproduction, carry, habitat, bevmort), .fns = list(mean = mean, sd = sd), na.rm = TRUE), .groups = 'drop') %>%
  dplyr::mutate(across(where(is.numeric), round, 3)) %>% 
  dplyr::filter(t>25) #remove burn-in

# write table into .csv file
write.csv(mean_stDEV, "mean_values_31Aug.csv", row.names=FALSE)
invisible(gc())

Abundance through time

abund_plot
## Don't know how to automatically pick scale for object of type <integer64>.
## Defaulting to continuous.

datatable(abund_quant, caption = "Mean abundance values for 2100 (t=115) and 2011 (t=25) and quantification of the abundance change between 2100 and 2011")

Reproduction through time

reprod_plot

Habitat suitability through time

habitat_plot

Carrying capacity through time

carry_cap_plot

Mortality through time

mort_plot

Abundance mismatch

As of 30th August 2024, the abundance mismatch was calculated only for landscape cells considered suitable (meaning it was calculated only for cells where the habitat suitability in time step 115 increased in relation the habitat suitability in time step 25). Additionally only time step 115 was plotted.**

#################################
## ABUNDANCE MISMATCH BOXPLOTS ##
#################################

abund_mismatch_df <- all_data[which(all_data$habitat_change >0), ]
invisible(gc())

abund_mismatch_df <- abund_mismatch_df[abund_mismatch_df$t==115,]
invisible(gc())
#summary(abund_mismatch_df)

# abundance mismatch
abund_mismatch_df <- abund_mismatch_df %>%
  group_by(species, scenario) %>% 
  mutate(abund_mismatch <- carry - abundance) %>%
  ungroup() %>% 
  rename(abund_mismatch = 15)
invisible(gc())

Abundance mismatch (in suitable cells for time step 90) boxplots coloured by the mean values tendencies, with outliers and the 10% and 90% quantiles

abund_mismatch_plot_GWF_withoutliers
## Don't know how to automatically pick scale for object of type <integer64>.
## Defaulting to continuous.

Abundance mismatch (in suitable cells for time step 90) boxplots coloured by the mean values tendencies, without outliers and the 10% and 90% quantiles

abund_mismatch_plot_GWF
## Don't know how to automatically pick scale for object of type <integer64>.
## Defaulting to continuous.

Step 4 - Cell-level plots (per species)

##################################################
# STEP 4 - CREATE CELL LEVEL PLOTS (per species) #
##################################################

# abundance change proportion
all_data <- all_data %>%
  group_by(species, scenario) %>% 
  mutate(abund_change_prop <- (abundance - abundance[t == 25])/abundance[t == 25]) %>%
  rename(abund_change_prop = 15)  %>% 
  ungroup()
invisible(gc())


# subset data for the last time step (t = 90)
abund_t115<-all_data[all_data$t == 115, c(2,3,9,10,4,15)]

#get the names of each scenario
scenarios <- unique(abund_t115$scenario)

# Define the RdBu palette with a midpoint at 0
palette <- brewer.pal(11, "RdBu")

These plots were created with for loops to allow for an image with a map per species for each scenario. In total there are 7 images with 9 maps each.

Plot for Scenario RCP2.6

### Plot for Scenario SSP1-RCP2.6

### Plot for Scenario SSP1

### Plot for Scenario RCP8.5

### Plot for Scenario SSP5-RCP8.5

### Plot for Scenario SSP5

### Plot for Scenario control

datatable(abund_chnage_prop_quant, caption = " Quantification of the proportion of abundance chnage between 2100 and 2011")

An extra plot was created representing Spatial hotspots of change in the landscape across all 9 species (MPCAS = Mean proportion of change across species).

spatial_hotspots_map

correlation_plot
## Warning: Removed 185535 rows containing non-finite outside the scale range
## (`stat_smooth()`).

Step 5 - Sensitivity analyses (per species)

Check this script to see how the sensitivity analysis was done.

See here image from 29th July 2024.